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We use quantum Monte Carlo simulations to study the phase diagram of hard-core bosons with 
short-ranged attractive interactions, in the presence of uniform diagonal disorder. It is shown that 
moderate disorder stabilizes a glassy superfluid phase in a range of values of the attractive interaction 
for which the system is a Mott insulator, in the absence of disorder. A transition to an insulating 
Bose glass phase occurs as the strength of the disorder or interactions increases. 



I. INTRODUCTION 

The interplay of superfluidity (SF) and localization in 
disordered Bose systems has been the subject of intense 
study for two decade^. Most of the theoretical investiga- 
tive effort has been focused on the quantum phase tran- 
sition between a superfluid and insulating phase. For 
example, it has been established that disorder leads to 
the appearance of an (insulating) Bose glass, sandwiched 
between the superfluid and Mott insulating phases. This 
topic is enjoying continued interest, especially since cold 
atom physicists have rec entl y produced controllable dis- 
order using laser speckled^ and looked at such phenom- 
ena as An derson localization of a one-dimensional con- 
densation! and the suppression of the condensate fraction 
in three dimensions-l^ 

The idea of "superglass" has come to the forefront in 
the context of the investigation of the (super)solid *He. 
The superglass phase was initially observed in quantum 
Monte Carlo simulations, in which the superfluid phase 
had an inhomogeneous condensate map on a microscopic 
scale.^ Biroli et al. proved that such a superglass phase 
does exist (at least as a metastable phase) by introducing 
an (artificial) model which could be mapped to a clas- 
sical system of hard spheres and studied in a controlled 
fashion.^ Recent experiments on solid *He have confirmed 
the strong interplay between a superfiuid component and 
a slow (glassy) dissipative component.'^ However, little is 
known yet about the superglass phase, and specifically 
about the actual role of disorder in promoting or en- 
hancing superfluidity. Give the current controversies and 
puzzles surrounding the interpretation of experiments on 
the possible supersolid phase of helium, further investi- 
gations of superfluid glassy phases are thus warranted. 

In this work, we provide strong numerical evidence for 
disorder-induced superfluidity in a lattice realization of 
hard-core bosons with a strong nearest-neighbor attrac- 
tion, in the presence of external disorder. In particular, 
we show that at low temperature and in a small range of 
attractive interactions, disorder of sufficient strength sta- 
bilizes a "glassy" superfluid phase. The superfluid den- 
sity reaches a maximum and then decays as the strength 
of the disorder increases, as an insulating glassy phase 
intervenes. In contrast to the case of repulsive bosons 
where disorder reduces the size of the superfluid phase. 



we see that strongly attracting hard-core bosons can be 
stabilized and made superfluid by disorder. In other 
words, disorder induces superflow in an otherwise insu- 
lating phase. Aside from supersolid ^He, such a scenario 
is possibly relevant to other condensed matter systems, 
e.g., high-temperature superconductors,'^ as well as to 
the elusive superfluid phase of molecular hydrogen,^ and 
to the role of substrate disorder in the superfluidity of 
(sub)monolayer helium films. 



II. MODEL 

We describe a disordered Bose system by means of the 
following Hamiltonian: 

H = -t'^{a\aj + h.c.) + V^hinj + ^5ihi . (1) 

i^j) (ij) i 

We consider here a square lattice oi N = L x L sites, 
with periodic boundary conditions. The sums (ij) run 
over all pairs of nearest-neighboring lattice sites, al (a.i) 
is the Bose creation (annihilation) operator for a particle 
at site i, hi — a^Oi is the local density operator. The first 
term of ([T|) represents the hopping of particles to nearest- 
neighboring sites. Henceforth, we choose the hopping 
integral t as our energy unit. 

The second term represents the interaction among 
bosons. A hard-core on-site repulsion is assumed, lim- 
iting the occupation of every site to no more than one 
particle. For the nearest-neighbor interaction, essentially 
all previous work based on ([T|) has focused on the repul- 
sive case, i.e. > 0, chiefly to elucidate the nature of 
the disorder-driven superfluid to insulator transition.!^ 
An enhancement of superfluidity by disorder has been 
predicted in some cases.fl^' Here, on the other hand, we 
consider the case of attractive nearest-neighbor interac- 
tion (i.e., negative V = —\V\), and neglect interactions 
among particles lying at distances greater than nearest 
neighbors. With this choice, the Hamiltonian (jlj, which 
is essentially a lattice model of quantum "sticky" spheres, 
is isomorphic to that of a spin- 1/2 XXZ quantum ferro- 
magnet. 

We model disorder by means of a random on-site po- 
tential Si, uniformly distributed in the interval [—A, A]. 
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Other theoretical representations of a disordered environ- 
ment could be considered, e.g., one in which the hopping 
matrix element t randomly varied from site to site, but 
in this work we restricted ourselves to the above, widely 
adopted diagonal model of disorder. In the spin lan- 
guage, the disordering potential is equivalent to a random 
on-site magnetic field along the z axis. 

In the absence of disorder, the ground state of ([T]) is 
a superfluid for \V\ < 2, whereas for \V\ > 2 only a 
Mott insulating phase exists, with exactly one particle 
per site, regardless of lattice geometry and dimension- 
ality. This is simply because the system can maximally 
lower its energy by having each particle surrounded by as 
many nearest-neighboring particles as possible, trumping 
any contribution from the hopping term. The regime of 
interest in this work is the latter, i.e., that in which no 
superfluid phase exists in the absence of disorder. 
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III. METHODOLOGY 

We perform grand-canonical quantum Monte Carlo 
simulations to study the ground state properties of (jl]), 
using the Worm Algorithm in the lattice path-integral 
representation.'i^^'* As the details of this computational 
method are extensively described elsewhere, and because 
the calculations performed here are standards, we shall 
not review it here, and simply refer interested readers to 
the original references. 

The results shown here correspond to a temperature T 
sufficiently low (typically /3 = 1/T = L), to be regarded 
as essentially ground state estimates. Simulations are 
carried out over square lattice of size varying from L=12 
to L—96, and estimates are averaged over a number M 
of independent realizations of the random disordering po- 
tential, typically M=100 (20) for L=12 (96). 



IV. RESULTS 

Fig. [l] shows the average particle density p as a func- 
tion of the chemical potential /i, for a particular value 
of |y| greater than 2 (|y|=2.3). For weak disorder (i.e., 
small A), the ground state of the system has exactly 
one particle per site, with an abrupt density jump at 
m/I^I = 2, when the lattice turns from empty to fully 
filled. However, for disorder of sufficient strength (figure 
shows results for A = 3), the density jump disappears, 
being replaced by a smooth curve, signaling continuous 
dependence of density on chemical potential. In other 
words, the disorder stabilizes phases at intermediate den- 
sities, consisting of interconnected "clusters" of particles, 
pinned by local fluctuations of the disordering potential. 
In this situation, the value = 2 corresponds to a 

particle density p = 0.5. 

Clearly, the issue immediately arises of whether such 
disordered phases may turn superfluid at low T, and 
what the nature would be of such a disordered superfluid 



Figure 1: (Color online). Ground state density p versus chem- 
ical potential p for |V^|=2.3, for weak (A=0.5, squares), inter- 
mediate (A=3.0, circles) and strong (A=9.0, triangles) dis- 
order. Results shown are for a lattice of size L—96, and are 
obtained by averaging over 20 independent realizations of the 
disorder. Statistical errors are smaller than symbol sizes. 



phase, simultaneously featuring broken translational in- 
variance. We investigated the occurrence of superfluid 
behavior by directly calculating the superfluid density pg 
(using the standard winding number estimator). Fig. [2] 
shows ps as a function of particle density p, in the limit 
0, for the one of the choices of model parameters 
of Fig. [l] namely A = 3 and \V\ = 2.3. The superfluid 
density increases from zero and reaches a maximum value 
at half filling, where approximately 12% of the system is 
superfluid. 

Obviously, numerical data such as those shown in Fig. 
|2]must be extrapolated to the L — > oo, in order for us to 
be able to male confidently the statement that superflu- 
idity observed in these systems is not merely a finite-size 
effect but survives in the thermodynamic limit. The in- 
set of Fig. [2] shows a typical extrapolation; estimates are 
shown for the superfluid density obtained for a fixed par- 
ticle density p = 0.5, on square lattices of different sizes 
(12, 24, 48 and 96), for |F|=2.3 and A=3. It is worth 
restating that these estimates are obtained by averag- 
ing results corresponding to several independent realiza- 
tions of the disordering potential. Based on results such 
as those shown in the inset of Fig. [2] we conclude that 
the superfluid signal remains finite in the thermodynamic 
limit. In general, we have observed that results obtained 
on a lattice with L=96 offer a close representation of the 
physics of the thermodynamic limit, at least in the range 
of parameters discussed here. 

Additional numerical evidence of superfluidity is ex- 
tracted from the behavior of the one-particle density ma- 
trix n{r), shown in Fig. [3] for the case \V\ = 2.3 and 
A = 3. The data are consistent with quasi long-range 
order, as expected in two dimensions, namely a power 
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Figure 2: (Color online) . Superfluid density ps versus particle 
density p for |1/|=2.3 and disorder strength A = 3. Statistical 
errors are smaller than symbol sizes. Results shown are for 
a square lattice with L = 96, and P=L, and are obtained 
by averaging over 20 independent realizations of the disorder. 
The solid line is a guide to the eye. Inset: Superfluid density 
for a fixed particle density p=0.5, computed on square lattices 
of varying size L. Extrapolation to infinite system size still 
gives a finite superfiuid density. 




Figure 3: (Color online). One-body density matrix n{r) for 
|V|=2.3 and disorder strength A = 3. Statistical errors are 
smaller than symbol sizes. Results shown are for a square 
lattice with L = 96, and P=L, and are obtained by averaging 
over 20 independent realizations of the disorder. Data show 
a weak power law decay of n(r) at long distances. 



law decay of n(r) at long distances. 

The observed superfluid phase is ostensibly induced by 
disorder, which stabilizes uniform phases of filling inter- 
mediate between zero and one. In order to gain further 
insight and in-depth understanding of the role of disor- 
der in actually promoting superfluidity, it is of interest to 
study the competition between the strength of disorder 
(A) and that of the attractive boson interaction {\V\). 



Figure 4: (Color online). Maximum value of the superfluid 
density ps (attained for p = 0.5) versus disorder strength A 
for different attractive interactions V. Statistical errors are 
smaller than symbol sizes. Solid lines are only meant to guide 
the eye. 



For definiteness, we consider the case of half flUing, cor- 
responding to a maximum in the superfluid density (for 
those systems for which superfluidity is observed). The 
same trends are also observed away from half filling. 

When the disorder is weak (A <C \V\), it cannot break 
apart clusters of particles, hence the system remains in- 
sulating, as shown in Fig. [4] or by the vanishing compress- 
ibility K = dp/dfi for A = 0.5 in Fig. [T| There are thus 
macroscopic domains (empty or fully filled) with hidden 
long-range order in the system.^^^E^ When the disorder 
becomes of the order of the attraction (A < \V\), sites 
and regions begin to appear throughout the system where 
the chemical potential is low enough to rip particles off 
the cluster, which then breaks down into large grains. 
These particles, however, are still largely localized in the 
vicinity of the energetically favorable sites created by dis- 
order, as the curve for A 2 in Fig. [5] shows. 

If we further increase the disorder strength, the grain 
size decreases to a microscopic scale, a (relatively) large 
fraction of the particles are delocalized, and superfluidity 
along interfaces (ridges) becomes possible, as also shown 
in previous numerical studies.^ This effect takes place 
essentially due to percolation. We can thus say that the 
disorder counters the insulating trend caused by the at- 
tractions, and actually makes the system superfluid. This 
is shown in Fig. |4] where we see a rather large super- 
fluid fraction as a function of disorder. Naturally, as 
the disorder strength is increased even further, insulat- 
ing glassy behavior re-appears, because the disorder is 
now so strong that it can block any superfluid path and 
localize particles, much as in the case of repulsive inter- 
actions. 

A similar scenario takes place on increasing the interac- 
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Figure 5: (Color online). Maximum value of the superfluid 
density ps{ attained for p=0.5) versus absolute value of at- 
tractive interaction \V\ for different diagonal disorder A at 
inverse temperature (3 — 96. Statistical errors are smaller 
than symbol sizes. Results shown are for a square lattice 
with L = 96. The solid lines are a guide to the eye. 

tion strength |y| at constant disorder bound A, as shown 
in Fig. [5] We have ah'eady explained the steep decay of 
the curve corresponding to A = 2 < |y| above, due to 
the lack of carriers. When A is greater than \V\, the 
disorder is sufficiently strong to destroy all macroscopic 
domains, and superfluidity can occur all over the sample. 
But as the disorder becomes stronger, it prevents parti- 
cle world lines from winding around the lattice. In the 
regime of strong disorder, both the disorder and the at- 
tractive interactions contribute to suppress superfluidity, 
as regions with nearly uniform chemical potential will be 
insulating due to the strong attraction, which pulls par- 
ticles together in such regions. 

It is worth noting that the above scenario is quite 
different from that of the repulsive disordered Bose- 
Hubbard model, where regions of uniform chemical po- 
tential are crucial for stabilizing locally a liquid phase, 
and thus the superfluid properties and the compressibil- 
ity of the whole system.^] We also note here that this 
insulating phase is compressible, as shown in Fig. [l] for 
A = 9t, which justifies the nomenclature "Bose glass" ^ 
The compressibility at half filling goes from zero in the 
phase segregated regime (no disorder) , to a large value in 
the superfluid phase, and decreases then monotonically 
over the Bose glass phase, when increasing the disorder 
bound at constant interaction strength. 

This "superglass" phase can be visualized through lo- 
cal superfluid density maps, shown in Fig. [6] for a partic- 
ular realization of disorder. The local value is obtained 
by statistically averaging local contributions to the total 
superfluid density (i.e., to the square of the winding num- 
ber), which, in the case shown in Fig. |6]for A = 3 and 
\V\ = 2.3, amounts to slightly less than 7%. We found 
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Figure 6: Map of the local superfluid density for a particular 
disorder realization, on a square lattice with L=96, \V\ = 2.3 
and A = 3. The total superfluid density ps equals 0.068(1) 
for this run. The white areas are small insulating grains, 
connected by superfluid interfaces. 

that the covariance between the superfluid density and 
the disordering potential is virtually zero, i.e., for these 
values of the parameters the physics is mostly driven by 
the attraction between bosons, consistent with the pic- 
ture given above in the case of strong disorder and strong 
attraction. 



V. CONCLUSIONS 

We have shown that disorder induces a superfluid 
phase in a lattice system of hard-core bosons with a 
strong nearest-neighbor attraction. While the system 
without disorder is an insulator of the ferromagnetic Ising 
type, the disorder can induce an inhomogeneous super- 
fluid (or superglass) phase (corresponding to in plane or- 
der in the spin parlance) in a window of interaction and 
disorder strengths. For stronger disorder bounds, the 
disorder and the attractive interactions work together to 
localize the particles. The coherence induced by disor- 
der might easily be observable in time-of-flight images 
for ultracold atoms or molecules. 

In the absence of disorder, the physics of our model 
is reminiscent of that of molecular para-hydrogen, long 
speculated to be a potential "second superfluid" , due to 
the light mass of its constituents (one half of that of he- 
lium atoms). On the other hand, superfluidity is not 
observed in para-hydrogen due to the strength of the in- 
termolecular potential, which causes the system to crys- 
talize at temperatures significantly above that at which 
Bose Condensation is expected to take place.I^ Recent 
numerical studied have shown that disorder ought not 
give rise to a superfluid phase of para-hydrogen. Based 
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on the results obtained in this work, we may argue that 
para- hydrogen may be a system too "deep" into the insu- 
lating regime (i.e., the effective value of \V\ is too large) 
for disorder to stabilize a superfluid phase. 

On the other hand, the results obtained here suggest 
that disorder may be responsible for the observation of 
superfluidity is helium films at coverages corresponding 
to less than a full monolayer. It should be noted, though, 
that helium films are expected to be superfiuid at "neg- 
ative" pre ssure, in fact all the way down to the spinodal 
density,^^] rendering them essentially different than what 
is discussed here. 
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